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Abstract. 

We reexamine the spin glass (SG) phase transition of the ± J Heisenberg models 
with and without the random anisotropy D in three dimensions (d = 3) using 
complementary two methods, i.e., (i) the defect energy method and (ii) the Monte 
Carlo method. We reveal that the conventional defect energy method is not convincing 
and propose a new method which considers the stiffness of the lattice itself. Using the 
method, we show that the stiffness exponent 9 has a positive value (9 > 0) even when 
D = 0. Considering the stiffness at finite temperatures, we obtain the SG phase 
transition temperature of Tsg ~ 0.19 J for D = 0. On the other hand, a large scale 
MC simulation shows that, in contrary to the previous results, a scaling plot of the SG 
susceptibility xsg for D = is obtained using almost the same transiton temperature of 
Tsg ~ 0.18J. Hence we believe that the SG phase transition occurs in the Heisenberg 
SG model in d = 3. 



PACS numbers: 75.50.Lk,02.70.Lq,05.50.+q 
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1. Introduction 

For a long time, it has been believed that the spin glass (SG) phase is realized in 
three dimensions (d = 3) for the Ising model PI 12] but not for the XY and Heisenberg 
models. [HI IH 03 El Hj Thus the SG phases observed in experiments were suggested to 
be realized due to anisotropy. [HJ E| However, numerical studies in the last decade have 
revealed that the SG phase might be more stable than what has been believed so far. 
In a long-range Ruderman-Kittel-Kasuya-Yoshida (RKKY) model, it was shown that 
the SG susceptibility exhibits a divergent singularity at a finite temperature, even when 
the anisotropy is absent .[THI ITT] This behavior has been attributed to the randomness of 
the spin position (site random model) rather than the long-range nature of the RKKY 
interaction. In fact, a short-range site random model composed of ferromagnetic spins 
and antiferromagnetic spins was shown to exhibit a long-range order phase characterized 
by the co-existence of a ferromagnetic and an antiferromagnetic orders. J2| On the other 
hand, for the XY and Heisenberg bond SG models, Kawamura and his coworkers took 
notice of chiralities of the spins and showed that a chiral glass (CG) phase transition 
occurs at a finite temperature Tqg 7^ 0, but the spin glass phase is still absent . [T3 ] IT U IT5] 
They insisted that an anisotropy mixes the chiral freedom and the spin freedom and 
the SG phase transition occurs at T SG (= T C q)- This view of the SG phase transition 
is quite attractive, because it gives a novel picture of the SG phase. That is, in the 
picture, the CG phase realizes in the real world, not the SG phase. However, their 
bases of the absence of the SG phase are obscure. Moreover, since the chirality is 
described by the spin variables, then the origin of the CG phase transition might be the 
usual SG phase transition. In fact, recent studies of the aging effects of the spin and 
the chirality autocorrelation functionsil6 i and the developments of the SG and the CG 
susceptibilities [T7] by means of a nonequilibrium relaxation method suggested that, if 
the CG phase transition occurs, the SG phase transition occurs at the same transition 
temperature T SG = T CG . Quite recently, Lee and Young presented the same conclusion 
using a finite size analysis of the correlation length of the spins and chiralities. |f 8j 

During the last decade, new algorithms for simulating the complex systems have 
developed and available computer power has increased enormously. It is therefore 
possible to reexamine in detail the SG phase transition of the Heisenberg model on the 
base of usual analyses. Here we consider Heisenberg models with and without random 
anisotropy on a simple cubic lattice described by 

H — — YliJijSiSj + £ DgS?Sf], (1) 

{ij) a^p 

where Si is the Heisenberg spin of \Si\ = 1 and Sf is its a-component (a = x, y, z), and 
(ij) runs over all nearest-neighbor pairs. The exchange interaction Jy takes on either 
+J or — J with the same probability of 1/2. We assume that the anisotropy comes from 
pseudo-dipolar couplings and impose the restriction D°f = Djf = D^ a . We further 
assume that D°f are uniform random values between — D and D. 

Evidences of the absence of the SG phase in the Heisenberg SG model which have 
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been believed so far are the following two points. 

(i) Negative stiffness exponent 9 at T = 0. [HI EJ [T3] 

(ii) Scaling plots of the SG susceptibility xsg an d absence of the crossing of the Binder 
ratio 0L-0IIS! 

Then we reexamine these two points to consider the possibility of the SG phase transition 
of the Heisenberg SG model. We will consider stiffness exponent 9 at T = and T/0 
in section 2, and properties of \sg and gi in section 3. We will give a special attention 
on the effect of the anisotropy, because it has been believed that the anisotropy brings 
the SG phase transition. So if it is true, we will find different properties between the 
models with and without the anisotropy. 

2. Stiffness exponent 9 

The most accepted evidence of the absence of the SG phase is results of the defect energy 
method. So we first consider the defect energy method. 

2.1. Conventional defect energy method 

The defect energy method comes from an application of a renormalization-group 
idea. [HI IU EOj That is, one evaluates an effective coupling Jl between block spins of the 
linear dimension L generated by the renormalization. The effective coupling Jl would 
depend on L as Jl ~ JL 6 with 9 being called as the stiffness exponent. When 9 > 0, 
the SG phase transition occurs at a finite temperature, while no phase transition occurs 
at any finite temperature when 9 < 0. To estimate Jl, one considers the domain wall 
energy AEl which is defined as the difference in the ground state energy of two lattices 
A and B of size L x L x L with the same bond distribution but with different boundary 
conditions. That is, for the lattice A, a periodic boundary condition is applied for every 
direction, and, for the lattice B, an antiperiodic boundary condition is applied for one 
direction and the periodic boundary condition for the other directions. By using the 
method, Banavar and Cieplak firstly estimated the value of 9 ~ — 1 and predicted that 
the SG phase transition occurs at T = 0. 3J Successive studies also predicted negative 

values for 9, i.e., 9 0.650 and 9 0.49. [Ej 

Recently, however, a doubt was thrown to the estimation of Jl.[2H[22] That is, in 
the calculation of AEl, one expects that no domain wall exists in the lattice A(or B) 
and hence one domain wall arises in the lattice B(or A). This expectation might be true, 
but another possibility would be equally true. That is, some domain wall will occur in 
the lattice A and some different domain wall in the lattice B. Then one might examine 
merely difference in the energy between those two domain walls. Does this defect energy 
AEl really gives the effective coupling Jl between the block spins? So we first examine 
AEl of the model with and without random anisotropy. We calculate AEl for lattices 
of L < 9, and for each L, the sample averages are taken over about 4000 independent 
bond realizations. Results of [| Ai£x|] are presented in Fig.l in a log-log form, where [• ■ ■] 
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Figure 1. The lattice size dependence of the naive defect energy [|A£^|] of the 
Heisenberg SG model in d = 3. 



means the sample average. Data for D = are curved. The most surprising thing is 
that this L dependence of [|A.Ex|] is similar to those in the case of D ^ 0. These results 
suggest two possibilities. One is that, a finite size effect is so large that the asymptotic 
region has not yet been reached. Since the curvature is upwards, it is possible that 
9 ~ or even more 9 > in the limit of L — > oo. The other is the inadequateness 
for estimating the defect energy as pointed out above. In order to examine the latter 
possibility, we study this problem in a different method. 

2.2. Stiffness of the system 

Apart from the renormalization-group concept, we consider the stability of the spin 
configuration of the system itself. The strategy of our examination is as 

follows. [2T| |2H I2SI HH3 We prepare a cubic lattice of L x L x (L + 1) with an open 
boundary condition in one direction of (L+ 1) lattice sites (z-direction) and the periodic 
boundary conditions in the other directions. That is, the lattice has two surfaces Q± 
and We call this system the reference system. First we determine the ground 

state spin configuration of the reference system. Hereafter, the ground state spin 
configurations on Qi and Ql+i are denoted as {S.^ 1 ^} and {S^ L+1 ^}, respectively. In 
this spin configuration, any distortion (domain-wall) in the z-direction will be removed, 
because the lattice has free surfaces Q\ and Ql+i- Then we add a distortion inside 
the system in the manner that {S^ 1 -*} are fixed and {S^ L+1 ^} are changed under the 
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condition that the relative angles between the spins are fixed. The ground state energy of 
this system is always higher than that of the reference system. When D = 0, this excess 
energy is the net one added inside the reference system, because the surface energy of 
JIl+i, which is given as the sum of the exchange energies between the spins on is 
conserved. We consider the stability of the system on the basis of this excess energy. One 
might think that the fixing of the relative spin directions on fij and Ql+i overestimates 
the stability of the spin configuration. We think, however, that this restriction is not 
serious for discussing the stability, because the increase of the excess energy to infinity 
for L — > oo means nothing but the existence of a strong correlation between the spin 
configurations on Qi and VIl+i- In fact, the same method was successfully applied to 
the Ising SG model in al = 2. gll 123 

We calculate two kinds of excess energies. One is the excess energy which is gained 
by rotating {S; (L+1) } by the same angle <fi around some common axis (z-axis) and the 
other is the excess energy which is gained by reversing {Si {L+1) }. Hereafter, we call the 
former system the rotated system and the latter system the reversed system. We think 
that it is sufficient to examine these two excess energies for considering the stiffness, 
because we can change {Si- L+1 ^} into any direction by combining the rotation and the 
reversal. The excess energy for the rotation AE rot ((j)) and that for the reversion AE Tev 
are given as 



where E G is the ground state energy of the reference system, and E TOt ((j)) and E rev are 
the ground state energies of the rotated system and reversed system, respectively. The 
lattice sizes studied here are L = 3 — 8 and, for each L, the sample averages are taken 
over about 1000 independent bond realizations. 

In Fig. 2, we present L- dependence of [AE rot (TT /2)]. Here we show data only for 
D = 0, because in the case of D ^ we could hardly evaluate the net excess energy of 
[AELt(7r/2)].[26j We clearly see that the data increase with L. From the slope of the 
asymptotic line shown in the figure, we tentatively determine the value of the stiffness 
exponent as 9 TOt ~ 0.62. That is, the SG phase would not be destroyed by a rotational 
perturbation. 

In Fig. 3, we present L-dependences of [A-E^v] for both D = and D ^ 0. Data 
depend little on the value D. For each of D's, they seem to lie on a curve with a common 
positive slope of 9 Tev ~ 0.4. Again, we get the view that the SG phase is stable at T ^ 
for both D ± and D = 0. 

Our results suggest 6 Tev , 9 rot > for both D = and D^0. However, their values 
are somewhat different with each other. It should be pointed out, however, that these 
values of 9 vev and 9 TOt may vary for L — ► oo, because in the lattice size range studied 
here [AE mt (ir/2)} is smaller than [A.E rev ] and the former increases more rapidly than the 
latter (0 rot > rev ), then as L increases further they would come close with each other. 
The convincing values of 9 rot and 9 rcv would be given in that range of L. Unfortunately, 




(2) 
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Figure 2. The lattice size dependence of the excess energy [A_E rot (7r/2)] of the 
Heisenberg SG model in d — 3. 



the lattice sizes are still small to examine this speculation. Any way, both the analyses 
of [AE TOt (ir / 2)] and [Ai? rev ] suggest that the system tends to be rigid as the size of 
the lattice becomes larger. Note that we have also calculated the defect energies of the 
system for D = using two replica boundary conditions [27] and found that they also 
increase with similar, positive slopes of ~ 0-59 and O^ffl ~ 0.46 for the ir rotation 
around the z-axis and the reversion, respectively. [2E] Hence, we conclude that the defect 
energy method never gives evidence of the phase transition at T SG = 0. 

What is the SG phase transition temperature Tgc? We may estimate Tsg by 
calculating the excess free energy [AT rot (T)] and [AT rev (T)][29j at finite temperatures. 
In fact, we have also calculated these quantities for D = 0.[2E] The result of [AF rev (T)] 
is shown in Fig. 4. It is seen that, at high temperatures [AF rev (T)] decreases with 
increasing L, whereas at low temperatures they increase with L. One estimates the 
phase transition temperature from the crossing temperature of the free energies for 
various lattice sizes L. In the present model, the crossing temperature Tl for the lattice 
sizes L and L + 1 shifts systematically to the low temperature side with increasing 
L. Then, we assumed that Tl decreases linearly with 1/L, and estimated Tl for 
L — ► oo as Tqo/ J — 0.188 ± 0.015. Note that the same extrapolation for [AF rot (T)] 
gave T^/ J = 0.192 ±0.015. Therefore we may conclude that, if the SG phase transition 
occurs, the transition temperature is T SG ~ 0.19 J. 




Figure 3. The lattice size dependence of the excess energy [AE Tev ] of the Heisenberg 
SG model in d = 3. 



3. Monte Carlo Simulation 

Now we reexamine the SG phase transition itself. Here we consider the model on a 
simple cubic lattice of L x L x (L+ 1)(= N) with skew boundary conditions along two L 
directions and a periodic boundary condition along the (L + 1) direction. We perform a 
MC simulation of the two-replica systems of the spins {Si} and {Tj} using an exchange 
MC algorithm [SO]- We calculate the order-parameter probability distribution Pl(?) of 

= [<%-«)>], (4) 

where (• • •) and [• • •] mean the thermal average and the bond distribution average, 
respectively. Here Q is the spin overlap defined by 



Q 



\ 



^E(rt 2 - (5) 

a,/3 



with q af3 = jjJ2iLi SfTf. Using Pl{q), we obtain two conventional SG quantities, i.e., 
the SG susceptibility xsg and the Binder parameter g(L,T) which are defined by 

Xsg = 3iV[(g 2 )], (6) 
^,T) = i(ll- 9 ML), (7) 




Figure 4. The lattice size dependence of the excess free energy [AF rev (T)] of the 
Heisenberg SG model for D = in d = 3. 28 Data at T — are [AE rcv ] presented in 
Fig.3 

where [{q n )} = J q n PL(q)dq. We examine the size and temperature dependences of these 
quantities both for D = and for D^O. The linear sizes of the lattice studied here are 
L = 5 ~ 19. Equilibration is checked by monitoring the stability of the results against 
at least two-times longer runs. The numbers of the samples are 480 for L = 5 ~ 9, 192 
for L = 11, 96 for L = 15, and 48 for L = 19. 

In Fig. 5, we show results of SG susceptibility xsg of the model without the 
anisotropy (D = 0). As the temperature is decreased, xsg f° r larger L increases rapidly. 
If the lower critical dimension di is less than the lattice dimension, di < 3, and the phase 
transition really occurs at T = T$q, the data for different L will be scaled as 

Xsg = L 2 -^F{L 1 I»{T-T sg )) 1 (8) 

where v is the exponent of the correlation length and rj is the exponent which describes 
the decay of the correlation function at T = Tsg- The scaling plots obtained by assuming 
Tsg 7^ and Tsg — are shown in Fig. 6. The scaling with Tsg 7^ works better than 
that with Tsg — 0, even if the data for the smallest size L = 5 are ignored in the 
latter. [HB] Note that in the previous scaling analysis [9J, Tsg = was estimated using the 
data for lattice of L = 7 — 15. Here, we use the data for a wider temperature range and 
add the data of the bigger lattice of L = 19. The phase transition temperature and the 
values of the critical exponent estimated here are Tsg/ J — 0.18 ± 0.01, v = 0.97 ± 0.05 
and r] = —0.1 ± 0.1. We should emphasize that this value of T S g is in good agreement 
with that estimated from the excess free energy of Tsg/ J ~ 0.19. It is noted, however, 
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Figure 5. Temperature dependences of the spin-glass susceptibility xsg of the ± J 
Heisenberg model in d — 3 for different sizes of the lattice at D = 0. 

that the possibility of T SG = is not ruled out from the scaling plot of Fig. 6(b), 
because in that case the temperature range of T > 0.2 J would be out of a critical 
region. [33] As the anisotropy is added, the transition temperature increases with D, i.e., 
T SG /J = 0.32 ± 0.03 for D = 0.2J, and T SG /J = 0.65 ± 0.05 for D = 1.0 J. 

The Binder parameter g(L,T) is the other quantity for examining the SG phase 
transition. It is believed that, if the SG phase transition occurs at T$q, g(L, T)'s for 
different L cross at Tsg- Contrary to our expectation, as shown in Fig. 7, they neither 
cross nor come together at T SG . This result seems to give the opposite view about the 
SG phase transition. However, the absence of the crossing of g(L, T) was also seen in the 
±J Heisenberg model in four dimensions (d = 4);3T] in which the SG phase transition 
is believed to occur at some finite temperature even when D — 0.0 If the absence of 
the crossing of g(L,T)'s for finite L says nothing about the SG phase transition, the 
same would be true when the anisotropy is present (D ^ 0). Then, we also calculate 
the Binder parameter of the model with D ^ 0. Here, since the system for D^0 has 
the inversion symmetry, we also consider the spin overlap of the diagonal components 
for which, in eq. (4), 

1 N 

Qdiag = "77 S i T i ( 9 ) 

iV 8=1 

is used instead of Q. Hereafter, we call the SG susceptibility and the Binder parameter 
calculated using Qdiag as the diagonal SG susceptibility and the diagonal Binder 
parameter and denote ^( diag ) and g^ dllLS \L,T), respectively. Results are presented in 
Figs. 8(a) and 8(b) for g(L,T) and g^ dmg \L,T), respectively. In fact, g(L,Tys for 
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Figure 6. Typical examples of the finite size scaling plot of the SG susceptibility at 
D = for (a) T SG ^ and (b) T SG = 0. 
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Figure 7. Temperature dependences of the Binder parameter g(L,T) of the ±J 
Heisenberg model for D = in d = 3 for different sizes of the lattice. The arrow 
indicates the transition temperature Tsg estimated from the scaling plot of xsg 

different L neither cross nor come together. In contrast, g(L,T) for larger L exhibits 
a dip. As D is increased, this property becomes more prominent. On the contrary, 
gr( dia s)(L, T) exhibits a usual behavior. That is, as the temperature is decreased, 
g( diag )(L, T) increases monotonically and its size dependence reverses. We suggest, 
hence, that the definition of the Binder parameter in terms of Qdiag is adequate for 
examining the phase transition for D^O and its crossing behavior supports the presence 
of the phase transition. It is noted, however, that the crossing temperature seems to 
deviate considerably from that estimated above. We think that this deviation comes 
from a finite size effect, because the crossing temperatures for different L's exhibit a 
considerable L dependence and, as L increases, it seems to approach T SG . We have also 
calculated ^ diag )(L, T) in the case of D = and found the absence of the crossing; 
behavior. We think this difference in the behavior of g( diag \L,T) comes from the 
occurrence of the drift of the whole system due to the 0(3) symmetry for D — 0. 
In fact, the diagonal SG susceptibility Xsg for = have been found to be much 
smaller than Xsg? whereas that for D ^ is larger than Xsg-IHH] We believe, hence, 
that the absence of the crossing of the usual Binder parameter g(L,T) will not say the 
absence of the phase transition of this system. We speculate that, even when D = 0, 
if the system becomes free from the drift, g( dmg \L, T) might exhibit a similar crossing 
behavior. 

Recently, it was proposed that the quantities A(L, T) and G(L, T) that measure the 
order-parameter fluctuations (OPF) exhibit the crossing behavior at Tsg even if g(L, T) 
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Figure 8. Temperature dependences of (a) the Binder parameter g(L, T) and (b) 
the diagonal Binder parameter </ dlag )(L,T) of the ±J Heisenberg model in d = 3 for 
different magnitude of the anisotropy D and for different sizes of the lattice. Open 
symbols are for D = 0.2J and closed ones for D = 1.0J. Arrows indicate the transition 
temperatures Tsg estimated from the scaling plot of xsg 
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does not: [321 EH! EH 

do) 

G(L T) _ go- 2 ) 2 } - m 2 (n) 

Then we also calculate A(L,T) and G(L,T) and examine their L-dependences. In Fig. 
9, we show G(L, T) for different L at D = and D = 0.2 J. When D = 0, T)'s for 
large L(> 9) seem to come together near Tsg- This property becomes more prominent 
in the anisotropic case of D = 0.2 J where the data for smaller L(= 5, 7) join. We have 
also seen that A(L, T) for both D = and D ^ show a somewhat different crossing 
behaviors at a temperature a little higher than Tsg- 



4. Conclusion 



We have reexamined the spin-glass (SG) phase transition of the ±J Heisenberg models 
with and without the random anisotropy D in three dimensions (d = 3). Attentions 
have been paid on the results of (i) the defect energy method and (ii) the Monte-Carlo 
method, because the evidences of the absence of the SG phase transition at a finite 
temperature have been given by these two methods. Our results have been summarized 
as follows. 

(i) The stiffness exponent 9: We have shown that the previous result of 9 < is not 
convincing because of the two reasons, i.e., (a) the meaning of the defect energy [| AEl\] 
in the conventional method is not clear, and (b) even if the method is meaningful, the 
plot of [|A£x|] as a function of L curves considerably. We have proposed a new method 
which considers the stiffness of the lattice itself. By using the method, we have shown 
that 9 > for D^O and the same is true for D = 0. Having considered the stiffness at 
finite temperatures, we have obtained the SG phase transition temperature T SG ~ 0.19 J 
for D = 0. 

(ii) The Monte Carlo method: A large scale simulation have enabled us to make a 
scaling plot of the SG susceptibility xsg which suggests the finite transiton temperature 
of T SG ~ 0.18 J for D = 0. The quantities G(L, T) and A(L, T) that measure the order- 
parameter fluctuations have exhibited a merging behavior near Tsg, but the Binder 
parameter g(L,T) has not exhibited the usual crossing behavior. However, analyses of 
the model with have suggested that the absence of the crossing of g(L, T)'s will 
not mean the absence of the SG order in this model. 

Our results have focused that, in contrary to the common belief, the SG phase 
transition occurs at a finite temperature. It should be noted that the two different 
method have given almost the the same SG transition temperature of T$q ~ 0.19 J. This 
value of the transition temperature is very close to that estimated from nonequilibrium 
properties of the model, i.e., Tsg ~ 0.19 J from the aging effect of the spin 
autocorrelation function |16j and Tsg ~ 0.21J from the nonequilibrium relaxation 
methodjEj. Quite recently, Lee and Young [TH] studied the Gaussian Heisenberg model, 
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Figure 9. Temperature dependences of the order-parameter fluctuations G(L,T) of 
the ± J Heisenberg model in d = 3 for different sizes of the lattice : (a) D = 0, and 
(b) D = 0.2J . The arrow indicates the transition temperature Tsg estimated from 
the scaling plot of xsg 
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and suggested T SG /J = 0.16 ± 0.02. This value of the transition temperature is also 
reasonably close to our value, considering the difference of the bond distribution. Hence 
we conclude that the model exhibits the SG phase transition even when the anisotropy 
is absent and its transition temperature is T SG ~ 0.19J. 

We give two comments. One might think that, for a larger D, the Ising values 
of the exponents {6 ~ 0.220J, and v ~ 2. and n = — 0.3[38j) should be recovered. 
However, we consider that this opinion is not necessarily true, because the random 
anisotropy in the model of eq.(l) is not uniaxial. Of course, we could not rule out the 
possibility that a finite-size effect masks true values. The other comment is that the 
strange behavior of g(L,T) of the Heisenberg SG model will come from the choice of 
the order parameter. When D ^ 0, the order parameter should be chosen as the sum 
of the diagonal components of the spin overlap, because only the inversion symmetry 
exists. We speculate that the same will be true in the isotropic case of D = 0, though 
the 0(3) symmetry recovers. To examine this speculation, we are currently making the 
simulation removing the uniform rotation of the system. JS] 

The authors would like to thank Dr. T. Nakamura for his valuable discussions. 
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